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We present a simple formalism for the calculation of the derivatives of the electronic density 
matrix p at any order, within density functional theory. Our approach, contrary to previous ones, 
is not based on the perturbative expansion of the Kohn-Sham wavefunctions. It has the following 
advantages: (i) it allows a simple derivation for the expression for the high order derivatives of p; (ii) 
in extended insulators, the treatment of uniform-electric-field perturbations and of the polarization 
derivatives is straightforward. 
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I. INTRODUCTION 

Linear response methodsfi*^ within the density func- 
tional theory approach (DFT)f^ have been successfully 
applied to compute a wide range of properties in real ma- 
terials such as phonon dispersions, dielectric constants, 
effective charges^ii^ and NMR spectrai^ 

Beyond linear response, perturbation theory applied 
to the Kohn-Sham (KS) orbitals allows the calculation 
of the derivatives of the energy at any order This kind 
of approach has two disadvantages: (i) although the fi- 
nal result is gauge invariant, i.e. invariant with respect 
to an arbitrary unitary rotation in the space of the occu- 
pied KS-orbitalSf^ the formulation of the theory depends 
on the chosen gauge. This becomes apparent in the ap- 
plication of the KS-orbitals orthonormality constraints 
at high order, (ii) In the case of periodic systems, the 
treatment of a perturbation due to a uniform electric 
field is not trivial, because the position operator, nec- 
essary to describe such a perturbation, is ill-defined in 
periodic boundary conditions. Much effort has been de- 
voted throughout the years to overcome this last prob- 
lem. Early treatments of the electric field perturbation 
for the calculation of the second and third order suscep- 
tibilities are particularly complex.^ A simpler formalism 
for the calculation of the second order susceptibility was 
obtained in Ref. 13, taking advantage of the 2n + 1 the- 
orem and of a Wannier representation of the orbitals. 
Only very recently Nufies and Gonza^ were able to give 
an expression for the derivatives of the DFT energy, with 
respect to uniform electric fields, at any order, by intro- 
ducing in the Hamiltonian an additional term depending 
on the polarization Berry phasCf^ 

We remark that, although a perturbation due to a 
macroscopic uniform electric field is ill-defined on an in- 
dividual Bloch states, such a perturbation is well defined 
on individual Wannier states^iiS which can be obtained 
by a different choice of gauge. This consideration sug- 
gest that the two problems, mentioned in the previous 
paragraph, are related, and that both problems might 
possibly disappear using a perturbative approach which 
is not based on the perturbative series of the single KS 
orbitals, but is solely based on the properties of the elec- 
tronic density matrix p, which is a gauge independent 



operator. 

In a recent paper— we gave an expression for the sec- 
ond order derivative of p which allowed the efficient com- 
putation of Raman spectral In the present paper we 
derive a general expression for the n-th order derivative 
of p, using the two relations = p, and [p, H] — 0, being 
[, ] a commutator and H the KS Hamiltonian. 

To fix our notation we define the electronic density 
matrix as 

p = X! \'^-")(^'"\' 

V 

where, throughout the paper, v or v' is an index running 
on the occupied valence states, and l^py) are normalized 
KS-eigenstates, i.e. H\ipy) = ey\ipy). Given a perturba- 
tion associated with a small parameter A, for a generic 
quantity F, we consider the perturbation series: 



F(A) = + AF(1' + A2f^(2) _^ ;^3^(3) _^ 



(1) 



The generalization to the case of different perturbations 
Ai, . . . A„ is straightforward, p and H stand for p(A) and 
H{\). We call Py, and Pc, respectively, the projec- 
tors on valence and conduction band states, i.e. Py = 



Given an Hermitean operator A we 
PvAPv,Acv = PcAPv, 



p(«),Pc = l-p(°^. 
define Acc = Pc^Pc, ^vv 
and Ave = {AcvV ■ 

The work is organized as follows. In Sec. |nl we use 
the relation = p to express p*^"' as a function of the 
operators p^y, having i < n. In Sec. lIIII we use the rela- 
tion [H, p] = to obtain an expression for p^y that can 
be easily computed using standard linear response tech- 
niques. In Sec. IIVI we show that, within our formalism, 
the perturbations due to a uniform electric field are well 
defined in extended insulators. In Sec. we derive a 
simple expression for the derivatives of the polarization. 



II. AS A FUNCTION OF {p^V}> WITH i < n 



We decompose p in pcc + Pvv + Pcv + Pvc, and we 
consider these four terms separately. The idempotency 
condition, p^ — p, implies that PcpPc = PcPpPc = 



2 



Pcp{Pc + Pv)pPc, or 

Pec — Pec Pec 



PcvPvc- 



When all the eigenvalues of pcc are lower than 1/2, i.e. 
for A sufficiently small, this relation between the two op- 
erators Pcc and pcvPvc can be inverted to obtain: 



Pcc 



1 - VI - '^PCVPVC 



(2) 



where the right-hand side denotes the operator obtained 
substituting pcvPvc in the Taylor series 



1 - VI -4a; 



2x^ + hx'^ 



(3) 



III. 



COMPUTATION OF 



In order to compute we introduce the wavefunc- 

IV'i*''') being an unperturbed 



tion |7?i"^) = Pcp(")|4°^) 
KS eigenvector. We have: 



(n) 

Pcv 



Y,Pcp^-^i^^:^){^^: 



(0)1 



E 



Equating to zero the n-th order term of the perturbation 
series of [H, p] = 0, we find: 



In a similar way, defining Ap w = pvv — P^'^'' ■, 

ApyyApyy + Apw = PW PVV - PVV = ^PVCPCV- 

When all the eigenvalues of Apyy are larger than —1/2, 
i.e. for A sufficiently small, lS.pvv can be expressed as a 
function of PvcPcv'- 



Apvy = Pvv 



1 - VI - 4pycpcv 



(4) 



Finally, p*^"^ can be expressed as a function of the 
{/^cy}' "with * < using the relation 

P^ ' =Pcv + Pvc + Pcc + Pvv^ (5) 

and taking the n-th order variation of Eq. (0) and Eq. Q 

through Eq. (j^J. As examples, observing that = 0, 
is easy to show that 



Note that each p^y is a gauge independent operator. 

In Eqs. 10-10, p*-"^ is expressed as p^^y + py^, plus a 
commutator, for n < 4. This property is used in Sec. |V] 
to compute the derivatives of the polarization. It holds 
at any order n. Indeed, as we show in the appendix: 



A) 


= f^l + 


P^l 




(6) 


.(2) 


= P'^ + 


Pfc^ 


' bcv' Pv^c\ 


(7) 


.(3) 


= P^^ + 


pfc^ 


' [Pcyj Pv^c\ + [p'cv^ Pv\j\ 


(8) 


.(4) 


- p'c% + 


P%^ 


' [Pcy' Pv^c\ + [Pcvi Pyc] " 


h 




+ [p^cUp%] 


+ [p\jvPv]jp\jVT Pv]j\- 


(9) 



^[77«,p("-*)] =0. 



i=0 



Multiplying this relation on the left by Pc and applying 
to IV'u"') to the right, we derive: 



n 



Solving the linear system of Eq. H13|l one can obtain jry^,"') 
and, thus, p*^"'. Since the right-hand side of Eq. \Yi\ 
depends on that in turn depends on p^"-*, the sys- 

tem is to be solved self-consistently, e.g. by using an 
iterative procedure. Eq. I|13(l needs to be solved only 

for a finite number of |77w"'') functions, running the in- 
dex V on the sole valence states. The linear system of 
Eq. H13I) is analogous to the one that is to be solved 
in the standard density functional perturbation theory 
(DFPT),i^ thus Eqs. (|10ll!^ll3|l give an efficient algo- 
rithm that can be easily implemented in available DFPT 
codes (as the PWSCFi^ or the ABINITi^ code), to com- 
pute the derivatives of p at any order. 

Alternatively, Eq. ifT^ can be written as 




(0)\ 



(14) 



4n) 



(") 

Pcv 



where n > 2, and 
Ovc = 



in) 

Pvc 



i=l 
n-1 



cv^^vc J 



P^SUp'vUY.PcV\p%]^ 



Pvc- 



1 - VI - '^PcvPvc 



'^PcvPvc 



1 - VI - "^Pvc pcv 
'^PvcPcv 



-pcv- 



(10) 



(11) 



where Gy = I]c lV'c°^)('0c°V(£^"'' ~ Cc"'') is the unper- 
turbed Green function operator projected on the con- 
duction band, and the sum J^c restricted to the empty 
conduction-band states. From Eq. p4l) one can recognize 

that \r]i^^) = Pcl'ipv'^) and that jr^f,^') is the projection on 
the conduction band of the KS orbital v in the parallel- 
transport gauge of Ref. However, at higher orders, 
there is not such a simple relation between the jryi*'') func- 
tions, defined in the present paper, and the variations of 
the KS-orbitals. 

Finally, as examples, we write p*-*' as a function of the 



(0) 



3 



'qif') wavefunctions, for the three lowest order: 



V 
V 

-T.\^i''){v'^^\vi'W:'\ (16) 

v,v' 

\ V V 

v,v' J 

We already used Eq. ((T?))l in Ref. O, to compute the 
Raman tensor. 



IV. TREATMENT OF THE ELECTRIC FIELDS 

Thanks to the commutators in Eq. (|13() . all the quan- 
tities needed to compute p^"^^ are well defined in an ex- 
tended insulator, even if the perturbation A is the com- 
ponent Ea of a uniform electric field, i.e., if H'^^^ — 
—era + dV^'^'^ I ^Eot^ being the a*'' Cartesian com- 
ponent of the position operator r, e the electron charge, 
and y^^'^ the Hatree and exchange-correlation potential. 
In particular, in an insulator, the commutator [r, p^""^^], 
which appears in Eq. H13|l . is a well-defined bounded 
operator, since the variation of the density matrix is 
localized ((r"|p("~-'^)|r') goes to zero exponentially for 
|r"-r'| ^oo). 

To prove the localized nature of p*-""^-* in a peri- 
odic system, we notice that p^""^) can be written (see 
Eq. pSfl ') as a sum of operators of the type 

where |aki;) ^nd l^^ki;) are Bloch wavefunction, i.e. 
Wv) = e^''--|ak.)/VA^ and |/?k.) = e^''-'-|/3k„)y ViV, be- 
ing N the number of unit cells, |aki;) and |/3kt;) wave- 
functions periodic in the lattice, normalized on the unit 
cell. In an insulator, the operators 



are analytic in k and periodic in the reciprocal space. 
Cloizeaux has shown in Ref. that an operator having 
the properties of D is exponentially localized. 

The representation of p^'^~^^ in terms of D is also useful 
to obtain a practical expression for the calculation of the 
[r, p'-"^^-'] commutator. In the limit of a converged k- 
points grid, 

k 

since the integral over its period of the derivative of a 
periodic analytic function is zero, fl^ is the unit-cell vol- 
ume. From this it can be easily demonstrated that 



kv 



. d\akv){Pkv\ _i 
dka 



(18) 



The terms required in Eq. when the perturbation is 
a miiform electric field, can thus be computed using 



J,(0)i^l«k."')(/3k^'||J.(0) 



dka 



where IV'k*^,^) are the periodic part of the Bloch wavefunc- 
tions, i.e. l-^k!)^) = e*'^'""|V'^°^)/A//V, and the bra-ket prod- 
ucts on the right-hand side are performed on the unit cell. 
In practical implementation, the derivative with respect 
to ka in the right-hand side of Eq. (|19|l can be computed 
numerically by finite-differentiation, using an expression 
independent from the arbitrary wavefunction-phase, as 
in Refs. HlR 



V. DERIVATIVES OF THE POLARIZATION 

Finally, with the present formalism, the computation 
of the n-th order variation of the polarization density 
P'"', becomes natural. The components of P*-"^ can be 
written as: 



p(n) 



2e 



rr{r„p(")} 



(20) 



akti)(/3k«| 



where the factor two accounts for the spin degeneracy 
and Tr{A] is the trace of the operator A. We substitute 
from Eq. (|10|l into Eq. (|20|l . Using the trace property 
Tr{[A, B]C} = Tr{A[B, C]}, is easy to arrive at: 



p(n) 



2e 
2e 



Tr \ K,p^">-- - [va^P^pl^i. + - E ([^-PcvlOyc^^ - K,P%]0. 



(n-i) 
CV 



i=l 
n-1 



(21) 



4 



where n > 2, Tm{z) is the imaginary part of the complex 
number z, and we have written the operators Oy|^ as 



k 



Pvv 



o(0) 



-OvcPc'v = —pvcOcv- For n > 2 the 



71-th order variation of pcc and pvv are: 

n—l n — 1 



being Ixb^O^i^l^. 



VI. CONCLUSIONS 



(») 



(i) ^(n-i) 



Concluding, we presented a formalism for the calcula- 
tion of the derivatives of the electronic density matrix at 
any order, within the density functional theory approach. 
Beside being simple, this formalism allows the treatment 
of extended systems in the presence of an external uni- 
form electric fields in a natural way, without introducing 
in the Hamiltonian an additional term depending on the 
polarization Berry-phase. 



APPENDIX 

The operators defined in Eq. are well defined for 
A sufficiently small, since the series 



1 - VI -4a; 
2x 



l + x + 2x'^ + 5x^ 



has the same convergence properties of Eq. From 
Eq. of the text, pcc = PcvOvc = OcvPvc and 



Eq. H1(J|) of the text easily follows 

Writii 
we have 



Writing O'yl as a function of , at the lowest orders 



*-^yc — Pvc 
'^vc — Pvc 



o 



(3) 

vc 



O 



vc 



^vc — 



(3) (1) (1) (1) 

Pvc PvcPcvPvc 



(4) _ „(4) 



Pv'c + E Pv'cP^cv Pv'^c^^+j+'kA + 

^(5) , (t) (]) (k) r 

Pvc + 2^ PvcPcvPvC°^+j+k,5 + 

'^PvcPcv PvcPcvPvc^ 



where ^ is a sum on positive integers. These equations 
allows to compute p*-"', with n <6. 
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